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1 Introduction 



Quantum billiards have been extensively studied as model systems having a chaotic or integrable 
classical dynamics (see e.g. [jij). They have also been considered as simple models for atomic nuclei or 
metal clusters. More recently they have been studied experimentally in ballistic microstructures (2], H] 
and in microwave cavities ||. In most of the experimental studies the shape of the equivalent billiard 
is determined only on the average. For instance, metal clusters have an underlying ionic background 
which implies an unavoidable degree of roughness of order of the interatomic distance, i. e. of order 
of the Fermi wavelength. For microstructures the roughness is mainly due to the irregularities in the 
depletion layer which can be estimated of being also of order of the Fermi wavelength []|]. Moreover 
clusters are produced in large amounts in molecular beams and one has to consider an ensemble of 
shapes. In the same line one also frequently considers an ensemble of microstructures (typically ~ 10 5 ) 
with a size dispersion ranging from 2% §] up to 30% §. 

The mean free path in the experiments quoted above is larger than the typical distances in the 
system, thus the billiard model is still meaningful. Nevertheless it should be corrected due to shape 
irregularities. In this paper we make an attempt of studying this phenomenon by considering an 
ensemble of billiards (in any dimension) obtained by a random deviation from a fixed initial shape 
(hereafter denoted as the perfect or unperturbed shape). We speak below of rough or bumpy billiards. 

The paper is organized as follows : in Sec. 2 and 3 we derive a semiclassical trace formula averaged 
over the ensemble of rough billiards. As an illustration the method is applied in Sec. 4 to study the 
supershell oscillations in rough metal clusters. We give our conclusions and compare with previous 
works in Sec. 5. 

2 Green function in presence of shape disorder 

Modern semi-classics have made extensive use of trace formulae such as derived by Gutzwiller in 
the context of quantum chaos (see [jjj and references therein) or for quantum billiards by Balian and 
Bloch . In this approaches the level density is obtained by computing the trace of the Green function 
G(r B ,r A ,k), solution of the Helmholtz equation with Dirichlet or Neumann boundary conditions. It 
is written in the semi-classical limit as a sum of contributions arising from the classical orbits of the 
system. In the case of a billiard it reads schematically : 

G(r B ,r A ,k) = J2 V(ky( kL -^/2) . (1) 
A->B 

where the sum is taken over the classical trajectories going from point A to point B. In (|l|) k is the 
wave- vector, it is related to the energy by fik = y/2mE (or fik = E/c for microwave cavities) ; {i is a 
Maslov index and T>(k) an amplitude characterizing the trajectory of length L considered. A general 
expression of V can be found in [7| . 

Let us now treat the case of the rough billiard. We consider that the shape disorder is weak enough 
so that a point r c on the frontier of the bumpy billiard can be written unambiguously as : 

r c = r Co + h(r Co ) n , (2) 
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where r Co is a point of the frontier of the perfect billiard and no the normal at this point, h is a 
random displacement, the characteristics of which will be specified later. 

Let us branch the perturbation (|2[) starting from a perfect billiard. The direct orbit going from 
A to B without bouncing on the boundary is not affected. If the shape modification is small enough 
(this will be made mathematically precise later) orbits experiencing only a few bounces will not be 
drastically altered (see Fig. 1). At first order in the semi-classical approximation one will consider 
that only their change in length is of importance, because it appears in the rapidly oscillating term 
exp(ikL) of (|lj). The modification of the slowly varying amplitude T>(k) is simply neglected. Long 
orbits on the contrary experience many bounces and they may be completely different in the rough 
enclosure and the perfect billiard. They also will be drastically different from a bumpy billiard to 
another, and ensemble averaging will very efficiently damp their contribution. Hence it is legitimate 
in the semi-classical limit to work in a perturbative approach where the only extra contribution with 
respect to the perfect billiard is the modification of the length of the orbits in ([l]). The spirit of the 
present approach is very common (see e.g. || || [l(J) and the results are similar to those obtained by 
other techniques (see below). 

If we denote by 5L the modification of length of an orbit due to surface roughness, using the 
notations of Fig. 1, for a single bounce trajectory one obtains : 

5L = (AE + EB) - (AC + C B) ~ (AC + CB) - (AC + C B) . (3) 

In (|D Co (resp. E) is the point of specular reflection in the perfect (resp. bumpy) billiard. Since 
AE + EB is an extremum of the length, at first order it can be computed replacing E by a nearby 
point. This has been done in the r.h.s. of (^) where point C was used, C being the intersection of the 
normal to the perfect billiard at Co with the bumpy frontier (see Fig. 1). This can be easily extended 
to orbits with n bounces and simple geometry yields 

n 

<5L~2 5Xr Oi )cos0 o . , (4) 

3=1 

where the sum is extended over all the bouncing points Cj of the classical trajectory on the boundary 
of the perfect billiard (Q c . is the normal angle of incidence at point Cj, see Fig. 1). 
Then the semi-classical Green function in the rough billiard is written as 

G(r B ,r A ,k) ~ ^ O P(^( fcI '-/ i7r /2)exp(2iA ; ^Mrc J )co S c ,) , (5) 
A^B \ j=l J 

where the upper index indicates that the sum is taken over the trajectories of the perfect billiard. 
Careful derivation puts the following limitations to the use of Eqs. (|) and (|) : 

(a) \h\/L <C |V/i| <C 1, or in other words : L should be greater than the typical distance between 
two bumps, itself being greater than the amplitude of shape disorder. These restrictions insure 
that replacing E by C in (||) is legitimate and that 5L <C L. 

(b) k\h\ 2 <C L insures that using the approximate length L + 5L in the semi-classical formula yields 
corrections which are indeed small compared to the main term (|5|). 
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(c) One should also make sure that diffractive corrections to the leading order semi-classics can be 
safely neglected. Hence, the typical distance between two bumps (|/i|/|V/i|) should be larger 
than the wave-length (1/k). If not, the amplitude T>(k) is significantly modified by the surface 
roughness. 

Keeping in mind the physical examples given in the introduction one sees that among the above 
restrictions only the ones involving Vh are not trivially satisfied. Indeed in the case of a large shape 
disorder the distance between two bumps is of order of the amplitude of a bump (then V/i is of order 
1) and also diffractive effects will have to be taken into account. Hence (||) is rigorously applicable 
only for small roughness (characterized by the restrictions (a), (b) and (c)). Nevertheless we will see in 
section 3 that in this limit the effects of the surface roughness are already very noticeable. Hence one 
can argue that when (|) is no longer valid the associated oscillations in the level density are already 
almost completely damped (see (|T2|)). 

We recall that (|5|) is only valid for short orbits. The contribution of long orbits in the bumpy 
billiard cannot be inferred from the motion in the perfect billiard. In order to have a formula rigorously 
applicable let us now damp the contribution of long orbits by performing an ensemble average of the 
Green function. The computation is straightforward and the average quantity reads 

<G(r B r A ,k)>^ J2° V ( k y {kL ~^ /2) ftx(2fc cos 9 C .) , (6) 

where x ls the characteristic function of the rough shape [|ll]]. It is the Fourier transform of the 
probability density of the frontier displacement h : 

X (s) = / P (hy sh dh . (?) 

J — oo 

In @) we have made the hypothesis that the bounces were separated by a distance larger than 
the correlation length of the shape disorder, i.e. p(h(r Cl ),..., h(r Cn )) = p(h(r Cl )) x ... x p{h(r Cn )). 
This restriction is not necessary but it simplifies the presentation. Following Ref. JlO| one could 
think of deformations strongly violating this assumption : this would be the case for instance of an 
unperturbed circle going to a rough ellipse. In this case it will always be possible to use the present 
formalism if considering the perfect ellipse as the unperturbed billiard. 

Most of the authors (see [11|) choose a gaussian distribution for h with standard deviation a : 



1 h 2 s 2 a 2 
P( h ) = — 7H= ex P(-77T) and X(s)=exp( — ) . (8) 

Note here that on the basis of Eq. (||) only one cannot check the validity of the restrictions (a), 
(b) and (c) above : they mainly concern the correlation length of the random function h and not only 
the characteristic function x which is our unique ingredient. The restrictions will be fulfilled if a is 
smaller than the correlation length, itself being smaller than typical distances in the billiard. 

The physics embodied in Eq. (||) can be simply interpreted by noticing that the wave pro- 
pagates as in a perfect billiard with, at each bounce, an extra damping factor x(2fccos6? .) (i.e. 
exp(—2k 2 a 2 cos 2 9 C ) in the gaussian model (jq)). The general form xi^k cos 9 C ) is commonly obtained 



in Kirchhoff theory of wave scattering from rough surfaces [11]. As anticipated this damps very effi- 



ciently the contribution of orbits experiencing many reflections. The quantity ka cos 6 Ci is known as 
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the Rayleigh parameter Jl^ ] and characterizes to which extend an incident wave is sensitive to surface 
roughness : as one would intuitively expect, the sensitivity is maximum for perpendicular incidence 
(cos0 Cj ~ 1) and short wave lengths {ka S> 1). 



3 Trace formula in the bumpy billiard 

The next step in the derivation of a trace formula is to compute the level density by taking the trace 
of the Green function : 

p(k) =-^Im I d D r G(r, r, k) == Im Tr G(k) . (9) 

IT J IT 

In (||) the integral extends over the interior of the billiard and D is the dimension of space. G(k) is 
the operator whose matrix elements in configuration space give G(r B ,r A , k). d s accounts for a possible 
spin degeneracy (in this case d s = 2, d s = 1 otherwise). 

It is customary to separate the contribution of the quasi zero length orbits to which the semi- 
classical approximation (|l]) does not apply. These orbits contribute to the smooth part p{k) of the 
level density through the "Weyl expansion" (see e.g. [Hfl). In three dimensions and for Dirichlet 
boundary conditions it reads 

., , / Vk 2 Sk \ 

where V is the volume of the billiard and S its surface area. In the typical case of a bumpy sphere 
of radius R with disorder of type the average < p(k) > is easily computed ; one obtains < 
V >= 47r(i? 3 /3 + a 2 R) and < S >= 4ir(R 2 + a 2 ), a is supposed to be small compared to R, thus 
surface roughness poorly affects the smooth part of the spectrum : for practical computations we 
will approximate < p(k) > by the value p(k) in the perfect billiard. We have given here a generic 3 
dimensional example but the same holds in any dimension. 

A quantity of primary interest is the oscillatory part p(k) of the level density. As shown by 
Gutzwiller |jj and Balian and Bloch Q, inserting the semi-classical Green function (]!]) in Eq. (||) and 
performing a non trivial stationary phase analysis leads to a trace formula which reads schematically 



p{k) = ^^(fe)sin(feL + i/7r/2) . (11) 
PO 

The sum (11) extends over all the classical periodic orbits (PO's) of the system. As in (|T|) A(k) is 
a slowly varying amplitude and v a Malsov index characteristic of the PO of length L considered (see 
[|l], 0). Some PO's may form continuous families, i.e. some orbits - forming a continuous set - may all 
have the same length and the same topology (such as the bouncing ball orbit in the stadium billiard or 
the PO's in integrable enclosures). Two orbits of the same family differ only by their bouncing points. 
Each family is represented by a single term in the summation (11) but its amplitude is enhanced by 
additional powers of k with respect to the contribution of an isolated orbit. This is related in |l4j] to 
local (possibly global) continuous symmetries of the hamiltonian. 
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If one wishes now to write a trace formula for a rough billiard, contrary to what happens for the 
Green function, the form (|i~l| ) of p(k) is different in the bumpy and in the perfect billiard. The reason 
is that some PO's may appear in continuous families and roughening destroys these families. This 
can be illustrated on the following two dimensional example : consider a rectangular billiard of which 
one edge -the upper one say- has been modified to adopt a sinusoidal shape. One of the important 
continuous families of PO's in the perfect rectangle -the vertical bouncing ball- is reduced in the 
bumpy rectangle to only a couple of orbits (those hitting the sinusoidal upper edge at points with 
horizontal tangent). Nevertheless the level densities of the two systems are certainly very similar if 



the edge deformation is small. This type of problem has been recently addressed in Ref. [10, [15|] and 
deserves a careful treatment. For isolated orbits it might in some cases be explained semiclassicaly 
by the introduction of complex PO's in the trace formula (see the discussion in and also fl6[| ). 
However we can bypass this kind of subtleties when averaging over disorder because it is legitimate 
to permute the trace and the average : < Im Tr G(k) >= Im Tr < G(k) >. Hence < p{k) > can 
simply be computed by inserting the average < G(r, r, k) > in the trace (g). Since < G(r, r, k) > is 
evaluated by using the orbits of the perfect billiard, the saddle point can be performed in the usual 
manner and yields the average oscillating part of the level density : 



< p(k) >c±Y, Ak) sm(kL + vtt/2) ]J X (2k cos 9 0j ) . (12) 

PO j=l 

The index in the summation indicates as before that all the quantities are evaluated in the 



unperturbed billiard. Hence n in (12) is the number of bounces of a PO in the perfect billiard, the 6's 
are the normal angle of incidence. 

Formula ( |i~2| ) is the most important result of the paper. It is valid for rough billiards in any 
dimension. It is interesting to note that when considering an integrable perfect billiard with an 
ergodic perturbation, the contribution of short orbits is correctly accounted for by Eq. (|j~2|), even 
without ensemble averaging. By ergodic, we mean that "any statistical average taken over many 



different parts of one shape realization is the same as an average over many realizations" [11]. In the 
case of integrable unperturbed billiards, all the orbits occur in families and the spatial integration (||) 
over a continuous family is - by the hypothesis of ergodicity - equivalent to ensemble averaging (this 
was certainly not the case in the simplified example above of sinusoidal deformation of a rectangle 
billiard because sine is not an ergodic function) . We recall that this is not correct for long orbits which 
may be very different in the bumpy and the perfect billiard, in this case ensemble average is necessary 
to damp the associated oscillations. 

To fix the ideas we apply ( |l2"l ) to a bumpy sphere of radius R. The total average level density 
<p>~ p+ < p > is plotted on Fig. 2 for two values of the surface roughness in the gaussian model 
(||) : a = 0.04 R and 0.06 R. The amplitude A(k) for each PO in the perfect sphere can be found 
in |?J. 75 PO's were included, with a maximum length of 26 times the radius. More precisely in the 
terminology of |?J the maximum values of the parameters are t = 5 and p = 15 (t being the winding 
number of an orbit around the center and p the number of bounces). We also indicate with black 
arrows the location of the first 15 eigenlevels in the unperturbed sphere. One sees on the figure that 
in the lower part of the spectrum, the wavelength being large, surface irregularities do not perturb the 
eigenstates much and there is still a strong bunching of levels. This shell effect gradually disappears 
for increasing energies (when the wavelength becomes comparable with the amplitude of the disorder). 
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4 Supershells in rough metal clusters 



In this section we illustrate the above results with the example of shell structure in metal clusters. 
Rough clusters will be described by a simple model first introduced in Ref. [17] : N electrons are 
moving independently in a bumpy sphere. The radius R of the perfect sphere scales with iV so that 
the mean electronic density is kept constant and equal to its bulk value : R = rgiV 1 / 3 , r s being the 
Wigner-Seitz radius of the bulk material. The standard deviation from this average shape is of order of 
atomic distances i.e. of order r s . Note that such irregularities are to be taken into account even if the 
cluster is "liquid-like" : the mean velocity of the ionic cores is always by several orders of magnitude 
smaller than the typical electronic Fermi velocity. Hence, as far as electronic motion is concerned, 
the ionic cores can be considered as frozen and this automatically implies a certain degree of surface 
roughness. 

In the unperturbed sphere, due to the high symmetry of the potential there is a strong bunching of 
levels leading to shell structure and magic numbers (see e.g. fl8|| ). This shell structure is modulated 
when the cluster size grows (it first disappears, then increases, etc.). This is denoted "supershell 
structure" and was first noticed in Ref. Q and seen experimentally by the Copenhagen group [jOj ]. 
We will here study this effect in a rough cluster. 

One of the most important observables when studying shell structure is the shell energy which 
is the oscillating part of the total energy. Shells are experimentally detected on a mass spectrum ; 
roughly speaking clusters with relatively smaller total energy (corresponding to minima of the shell 
energy) are more stable and will be more numerous in a beam. This is also correlated with larger 
ionization potential, but for this observable shell effects decrease with cluster size, making its study 
more difficult. 

The Fermi wave vector kf and total energy Etot are defined by : 



N = M{k f ) 



p{k)dk and E tot (N) 



k f h 2 k 2 



2m 



p(k)dk , 



(13) 



io Jo 
In ( |l3| ) J\f{k) is the integrated level density, or spectral staircase. As the level density, J\f can be 
written as the sum of a smooth quantity N and an oscillating part J\f . The same holds for kf and 
E to t considered as functions of N. The smooth terms are defined by : 



N = /V{k f ) 



p(k)dk and E tot (N) 



h h 2 k 2 



p(k)dk . 



(14) 



/o Jo 2m 

As explained before we will identify the average smooth quantities < Etot > and < kf > with their 



value in the perfect sphere. On the basis of the Weyl expansion in the sphere and of Jl4| ) one obtains 
the following relations : 



kf(N) 



1 



Q2 
3a 3 



and 



ATX -1/3 
CL3 



3ai<23 



9ai 



N \ - 2 / 3 
d3 



(15) 



Et t(N) = e f 



303 

5 



02 (^\ 2I \<^_ 

2 \aj 



2(21(23 



3a 3 



ATX 1/3 
C*3 



(16) 
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where 03 = 2d s /97r, a2 = —d s /A, a\ = 2d s /3ir (see [|13|). Kf and Ef are the bulk Fermi wave-vector 
and Fermi energy : Kjr s = a 3 l ^ and Ef = h 2 K 2 /2m. 

The average oscillating part is obtained by subtracting (|l~6|) to the average total term. It is 
computed using the following approximation : 

~ f k f h 2 k 2 f <k f> h 2 k 2 

< E tot >=< / - — p{k)dk > -Etot ^ / - — < p(k) >dk- E tot . (17) 
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The r.h.s. of (0) can be considered as a simple first approximation of the exact result. It is made 
necessary by the difficulty mentioned above for estimating p(k) in an individual rough cluster. We 
can invoke the hypothesis of ergodicity of surface disorder to make this approximation sound : for a 
cluster such as created in a molecular beam one can argue that the individual level density will have 
a pattern very similar to the one displayed in Fig. 2 (c/the discussion in Sec. 3). Hence, although in 
a given cluster p(k) is exactly a sum of delta pics, the bunching of levels in an individual spectrum 
will disappear at the same wavelength that predicted in the average < p(k) >. 

Then < Etot > can be computed with the following scheme : one first determines < kf > by 
numerical inversion of the first term of Eq. ( |l3|) (with < p > replacing p) and the integral of the 
second term is then computed numerically. One can also evaluate the asymptotic behaviour of the 
integrals (^) to get an analytical expression. For this purpose the exact Fermi wave-vector in written 



as kf = kf + kf, where kf is given by (15) and kf is supposed to be small. Then performing on 



J\f = M + M a first order limited expansion one gets 

N{k f ) ~ N{if) +h{ d §)_ + + h _ • (is) 

\ / kf V / kt 

Since N(kf) = M(kf) = N by definition (see (|l3| , p^ )) one obtains : 

k t ~ - EM (19) 

This expression can be evaluated by keeping only the leading order of M (obtained by integration 
of (pT|)) and neglecting the oscillatory term at the denominator. Then averaging over disorder yields : 

< kf >^ -^rr E° co < k f L + ^/ 2 ) fl *(2fc/ c «s &j ) • (20) 

P\ K f) PO L j=l 

Then the shell energy Etot is computed from the difference between the total energy and its smooth 
part : 



~ f k f h 2 k 2 f k f h 2 k 2 

Etot = E tot ~ E tot = / ——p{k)dk + / -—p{k)dk 

Jk f 2m Jo 2m 

~ h 2 k 2 _ A n 2 k 2 ~ n 2 k 2 _ 

- "f-^m+l ^mdk +kf ^p(k f ). (2D 
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In (^) we replaced the integrals by their large N approximations. The dominant contribution is 
obtained by integrating by parts the integral appearing in the last term of (21) ; it cancels due to 
relation (|l9|). After averaging over disorder the next order reads : 



< E tot >~ < T(k f ) >~ ^ Y?~PT s[n C k f L + U7T / 2 ) f[x(2k f ^se Cj ) . (22) 

m lm pQ kfL j=1 

J-(k) is a primitive of J\f(k) which has been estimated at first order in the r.h.s. of Eq. (^2|). 
This formula gives an accurate approximation of the value of < Etot > computed numerically : the 
result is shown on Fig. 3 for the gaussian model <M) and for several values of the surface roughness 



(cr/r s = 0, 0.1, 0.2 and 0.3). The dashed line corresponds to Eq. (22) and the solid line to numerical 
evaluation of the integrals (|l3|]l7|) . Here the integration adds an extra smoothing compared to Fig. 2 
and one needs to take into account much less orbits, the figure includes 35 PO's up to a length L = 12R. 
Actually a very reasonable result can be obtained with only the 7 shortest orbits, we included more 
orbits here to have an accurate description of shell effect in the supershell region iV 1 / 3 ~ 8. 

We see on Fig. 3 that shell structure is very sensitive to surface irregularities of small amplitude. 
Nevertheless roughness reduces all oscillations without modifying the qualitative features of the su- 
pershells. Hence the present approach legitimates the usual explanation of supershell effects in rough 
metal clusters as being due to the interference of the square and triangular orbits 0, [l{J, although 
these orbits might not exist in an individual cluster. Including temperature effects and quantitatively 
comparing with the experimental results could fix an order of magnitude for the irregularities of the 
surface of large metal clusters. A very rough estimate based on separation energies for small clusters 
gives the value a ~ 0.2 r s 



5 Discussion 

In this paper we have studied the oscillating part of the level density in billiards with small size 
shape irregularities and we derived a semiclassical trace formula averaged over shape disorder. The 
important feature of the level density is the gradual disappearance of the oscillations with increasing 
energy : when the wavelength is of order of the typical size of the surface defects this induced a shift 
of the eigen-energies which leads after averaging to a structureless level density. 

The same type of approach has recently been used to study, within semiclassical approximation, the 
role of families of orbits broken by a small perturbation of an initial shape having a local continuous 
symmetry. The authors of Ref. [j31| study the bouncing ball orbit in a deformed stadium billiard 
and in Ref. [1C] a general trace formula is derived accounting for the role of broken families. In the 



present work, the main difference is the inclusion of an ensemble average yielding a formula valid also 
for isolated orbits. Moreover averaging damps the contribution of long orbits to which the simple 
perturbation technique (HJH) does not apply. This averaging method is motivated by the experimental 
techniques of mesoscopic and cluster physics. 

The method has been applied to study supershell oscillations in metal clusters using a model 
accounting for the irregularities of the surface due to the underlying ionic structure. It is also of 
interest in ballistics microstructures with shape irregularities [pif| . 
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The problem of surface roughness of metal clusters has recently been addressed using an other 
approach : in Ref. [22] the effect of disorder is represented via addition to the hamiltonian of a random 
matrix perturbation (c/Ref. |23|]). The results for the average level density and shell oscillations are 



qualitatively very similar to what is presented here. In addition the authors of Ref. 22] argue that 
these effects could explain experimental shifts in the measured mass distribution. Note that in Ref. p^] 
and also in the present study the effects of temperature are indirect : although the usual temperatures 
reached in experiments are small compared to the Fermi energy (one remains in the very degenerate 
limit k B T <C £f) they are sufficient to induce a disorder of the ionic arrangement which has a sizeable 
effect on shell structure. 
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Figure captions 



Figure 1. Classical orbits going from point A to point B with one bounce one the boundary. 
The thick lines correspond to the perfect billiard and the thin lines to the bumpy billiard. Point Cq 
(resp. point E) is the point of specular reflection on the perfect (resp. bumpy) enclosure. 6 Co is the 
normal reflection angle at Cq. C is the intersection of the normal at Cq with the bumpy boundary. 



Figure 2. Total average level density in bumpy sphere for a = 0.04 R (solid line) and a = 0.06 R 
(thick solid line). The dashed line represents the smooth term p(k). The black arrows indicate the 
location of the first 15 eigenlevels in the perfect sphere. 



Figure 3. Average oscillating part of the total energy is rough metal clusters as a function of 
N 1 / 3 for several values of surface roughness : a/r s = (upper graph), 0.1, 0.2 and 0.3 (lower graph). 
< E tot > is expressed in units of the bulk Fermi energy ej. The solid line corresponds to numerical 
evaluation of the integrals (|l^,17) and the dashed line to Eq. (22). 
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